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Research  Objectives 

To  develop  and  study  solution  algorithms  for  equations  used  to  model  reflection 
properties  of  underwater  objects  interacting  with  sonar  signals. 

Technical  Approach  and  Accomplishments 

The  emphasis  has  been  on  development  of  algorithms  for  the  numerical  solution 
of  the  Helmholtz  equation 

—  Au  —  k2u  =  /,  (1) 

the  solution  to  which  is  used  to  determine  the  acoustic  pressure  associated  with 
an  underwater  signal.  Here  we  consider  a  three-dimensional  box-shaped  domain 
fi  =  (ai,&i)  x  (02,^2)  x  (03,^3)  C  3f3,  with  Sommerfeld-like  boundary  conditions 

un  —  iku  =  0  (2) 

on  dtt. 

Discretization  of  the  problem  (l)-(2)  results  in  a  linear  system  of  equations 

Au  =  f.  (3) 

Since  the  problem  is  fully  three-dimensional,  any  reasonable  discretization  will 
contain  a  large  number  of  unknowns  and  require  considerable  storage.  Direct 
methods  based  on  Gaussian  elimination  with  partial  pivoting  require  a  prohibitive 
amount  of  additional  storage  and  thus  have  limited  use.  Multilevel  methods  suf¬ 
fer  from  the  requirement  that  the  coarsest  spaces  used  must  be  fine  enough  to 
accurately  represent  the  solution.  In  addition,  the  complex  symmetric  coefficient 
matrix  A  typically  has  eigenvalues  with  both  positive  and  negative  real  parts. 
This  can  cause  difficulties  for  iterative  solution  methods,  and  preconditioning  of 
the  matrix  is  essential  in  order  to  attain  efficiency. 

In  this  project,  we  solved  the  discrete  Helmholtz  equation  using  Krylov  sub¬ 
space  iterative  methods  with  a  preconditioning  methodology  derived  from  fast 
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direct  methods.  The  basic  principle  behind  fast  direct  solvers  is  to  apply  an  in¬ 
expensive  transformation  to  break  a  problem  into  a  number  of  lower-dimensional 
but  independent  problems.  Many  solvers  use  fast  Fourier  transforms  (FFT’s)  to 
achieve  separation  of  variables  and  then  solve  the  resulting  set  of  decoupled  prob¬ 
lems  using  sparse  matrix  methods.  Fast  direct  methods  are  standard  tools  for 
solving  the  Poisson  equation  on  regular  domains  with  Dirichlet,  Neumann  or  peri¬ 
odic  boundary  conditions;  they  can  be  adapted  to  other  domains  via  capacitance 
matrix  or  embedding  methods.  They  have  been  used  for  the  three-dimensional 
Helmholtz  equation  with  Dirichlet  or  Neumann  boundary  conditions  on  an  irreg¬ 
ular  domain,  and  for  the  two-dimensional  problem  in  polar  coordinates  with  non¬ 
reflecting  boundary  conditions  (derived  from  a  Dirichlet-to-Neumann  mapping). 
Here,  we  developed  efficient  solvers  for  problems  with  Sommerfeld-like  boundary 
conditions  on  box-shaped  domains.  Combining  our  techniques  with  capacitance 
matrix  methods  would  produce  solvers  for  general  geometries  in  Cartesian  coor¬ 
dinates,  including  exterior  problems. 

Our  accomplishments  were  threefold: 

1.  We  approximated  the  discrete  operator  A  with  a  matrix  Q  that  can  be 
treated  with  fast  direct  methods.  For  finite  difference  discretizations,  we 
derived  Q  by  defining  and  discretizing  the  differential  operator  in  the  same 
way  as  for  .4  except  that  the  boundary  conditions  on  either  two  or  four 
faces  of  Cl  are  replaced  by  more  convenient  ones  (Dirichlet  or  Neumann).  The 
resulting  matrix  Q  differs  from  A  by  a  (relatively)  low-rank  operator  and  can 
be  used  as  a  preconditioner  for  A,  to  accelerate  the  convergence  of  iterative 
solvers  based  on  Krylov  subspaces.  We  also  developed  variants  of  these  ideas 
for  finite  element  discretizations  (on  uniform  grids),  focusing  on  trilinear 
elements.  Here,  rather  than  explicitly  modifying  the  boundary  conditions 
to  construct  Q ,  we  used  the  fact  that  the  discrete  operator  A  is  close  to 
a  block  Toeplitz  matrix  and  replaced  certain  sub-blocks  of  A  by  Toeplitz 
approximations  that  are  amenable  to  fast  transforms.  For  both  types  of 
discretizations,  we  demonstrated  empirically  that  Q  meets  the  requirements 
for  an  effective  preconditioner: 

•  Applying  the  action  of  Q_1  to  a  vector  is  not  too  expensive.  For  our 
preconditioners,  using  Q~l  entails  a  set  of  FFT’s  together  with  solution 
of  smaller  dimensional  problems. 

•  Q  greatly  reduces  the  number  of  iterations  needed  by  Krylov  subspace 
methods  to  solve  (3). 

In  particular,  we  showed  that  for  several  choices  of  Q ,  the  experimental  con¬ 
vergence  behavior  of  preconditioned  restarted  GMRES  depends  onlv  mildly 
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on  both  the  wave  number  k  and  the  discretization  mesh  size.  In  addition,  we 
demonstrated  how  the  methods  can  be  implemented  on  a  parallel  computer 
with  high  efficiency. 

2.  We  then  worked  on  analyzing  this  algorithm  to  confirm  the  empirical  proper¬ 
ties.  We  explicitly  calculated  the  eigenvalues  of  the  preconditioned  operator. 
The  main  innovation  is  that  the  eigenvalues  for  two  and  three-dimensional 
domains  can  be  calculated  exactly  by  solving  a  set  of  one-dimensional  eigen¬ 
value  problems.  This  permits  analysis  of  quite  large  problems.  For  grids 
fine  enough  to  resolve  the  solution  for  a  given  wave  number,  preconditioning 
using  Neumann  boundary  conditions  yields  eigenvalues  that  are  uniformly 
bounded,  located  in  the  first  quadrant,  and  outside  the  unit  circle.  In  con¬ 
trast,  Dirichlet  boundary  conditions  yield  eigenvalues  that  approach  zero 
as  the  product  of  wave  number  with  the  mesh  size  is  decreased.  These 
eigenvalue  properties  yielded  the  first  insight  into  the  behavior  of  iterative 
methods  such  as  GMRES  applied  to  these  preconditioned  problems. 

3.  We  have  taken  our  algorithm  for  solving  the  discrete  indefinite  Helmholtz 
equation  on  a  three-dimensional  box-shaped  domain  with  Sommerfeld-like 
boundary  conditions  and  modified  it  to  make  it  efficient  for  inhomogeneous 
media  in  which  the  speed  of  wave  propagation  is  different  on  an  interior 
domain.  The  preconditioners  are  of  two  types.  The  first  is  a  precondi¬ 
tioner  that  we  designed  for  a  homogeneous  medium.  This  works  well  if  the 
volume  of  the  interior  domain  is  relatively  small  compared  to  the  overall 
volume  and  if  the  speeds  are  close  in  value.  The  second  preconditioner  uses 
a  direct  solver  on  the  interior  domain  in  addition  to  the  homogeneous  pre¬ 
conditioner.  We  presented  experimental  results  demonstrating  that  these 
algorithms  display  efficiency  comparable  to  that  we  showed  earlier  for  a  ho¬ 
mogeneous  medium,  making  them  quite  useful  for  the  problem  of  acoustic 
analysis  with  a  submarine. 


Relevance  to  the  Navy 

The  equations  of  structural  acoustics,  of  which  (1)  is  one  component  (to  be  coupled 
with  an  equation  of  elasticity  for  a  structure)  model  fluid-structure  interaction  in 
the  presence  of  an  acoustic  wave.  This  model  is  used  to  identify  the  location  and 
structural  properties  of  underwater  objects. 
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Table  1:  Representative  iteration  counts  for  GMRES(20)  with  preconditioners 
that  use  two  sets  of  one-dimensional  transforms,  for  finite  difference  discretization 
with  several  wave  numbers  and  grid  sizes. 
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Table  2:  Representative  CPU  times  for  solution  of  finite  element  discretizations, 
with  k  =  5  and  several  grid  sizes. 
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Representative  Results 

Our  results  are  in  the  form  of  iteration  counts  and  CPU  times  for  solving  the 
discrete  Helmholtz  equation  (1).  Representative  statistics  are  given  in  Table  1 
and  2.  Table  1  shows  iteration  counts  for  solving  the  discrete  problem  on  a  three- 
dimensional  m  x  m  x  m  grid  for  two  variants  of  the  preconditioner.  Numbers 
above  the  jagged  line  correspond  to  accurate  models  (with  at  least  ten  grid  points 
per  wave).  The  results  demonstrate  the  insensitivity  of  performance  to  mesh  size 
(1/m)  as  well  as  the  relatively  small  dependence  on  wave  number.  Table  2  shows 
CPU  times  on  an  IBM  SP-2  computer.  These  results  demonstrate  the  scalability 
of  the  algorithms.  We  note  that  the  largest  problems  that  we  have  solved  contain 
729,000  degrees  of  freedom  (on  a  90  x  90  x  90  grid)  and  required  just  two  minutes 
of  CPU  time. 
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